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Abstract 

The volume-of-tube formula was first introduced by ffotcUing (i939), to 
solve significance of terms in nonlinear regression models. Since this pioneer- 
ing paper, there has been significant work on extending the tube formula to 
more general settings, including multidimensional problems, and many new 
applications in statistical inference, including confidence bands in regression 
and smoothing models; applications to functional data analysis; testing in 
mixture models; and spatial scan analysis. 

Implementation of the tube formula requires numerical evaluation of cer- 
tain problem-specific geometric constants that appear in Hotelling's formula 
and its extensions. The purpose of this note is to describe a software library, 
libtube, that performs the calculations. A variety of illustrative examples 
are given. 

Source code for the libtube library and examples can be downloaded from 
http : / /www . herine . net/ stat/libtube/ . 

1 Introduction 

The volume-of-tube problem can be stated rather simply. Given a curve (or 
manifold) Ml lying in n-dimensional Euclidean space, what is the volume 
of the set of all points lying within a radius r of the curve? In statistical 
applications, the spherical version of this problem often arises; the manifold 
lies on the surface of the unit sphere in n dimensions, and one wishes to 
compute the (n — l)-dimensional volume (or surface area) of the set of points 
lying within a distance r of the manifold. 

The volume-of-tubes formula was formulated and solved by Hotclling (19-39), 
motivated by application to significance testing in nonparametric regression. 
A companion paper, Wcyl (1939), extended the results to higher dimensional 
manifolds; that is, when is a surface, or more generally when A4 is a, 
manifold of dimension d < n. 

The main purpose of this article is to describe a set of routines written by 
the author to implement the volume-of-tube formula in statistical problems. 
In section 2 the tube formula (with boundary corrections) is described. The 
libtube software is described in Section 3. Applications to non-linear regres- 
sion, simultaneous confidence bands and mixture modeling are described in 
Sections 4, 5 and 6 respectively. 
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Figure 1: The manifold is represented by the red curve. The tubular neighborhood 
is approximated by trapezoids, plus the two end-point caps. 

2 The Volume-of- Tubes Formula 

The volume-of-tubes formula was first derived by HotcUing (1939). The result 
can be illustrated on the plane by reference to Figure 1. The manifold is 
represented by the red curve. The tubular neighborhood of a given radius r 
is approximated by trapezoids, plus the two end-point caps. Adding up the 
area of the trapezoids and letting the partition become increasingly fine shows 
that the area (or two-dimensional volume) of the tube is 

Length of Manifold x 2r + Tvr^ . 

The 2r represents the cross-sectional area of the manifold, while ivr^ rep- 
resents the area of the end-point caps. The result extends to manifolds em- 
bedded in n-dimensional space; 

Volume = KoV„-ir"~^ + ^-^V„r". 

Here kq is the length of the manifold and is the number of end-points (often, 
lo = 2). The functions ilioir) are the cross-sectional area and volume of the 
end-point caps respectively, and 14 = n''^^ /T{1 + k/2) is the volume of the 
fc-dimensional unit sphere. 

When the manifold lies on the unit sphere, the result is similar, but the 
cross-sectional area is replaced by a certain partial beta function. The result 
is 

Volume = ■^^|^-P(-Bl,(„_2)/2 > w'^) + ^-^^^P{Bi/2.(n~l)/2 > li?), 

where Ba,b denotes a random variable following a beta distribution with pa- 
rameters a and h\ An — 27r"/'^/r(n/2) is the surface area of the unit sphere 
in 7?", and w = 1 - r'^ /2. 

Multidimensional Manifolds 

Figure 2 shows a tube around a two-dimensional manifold. To compute the 
volume of the tubular neighborhood, one divides the tube into different pieces: 
a main piece, the half-cylinders around each edge, and wedges at each corner 
of the manifold. For higher dimensional manifolds, the ideas are similar, but 
there are more pieces to take care of. 



2 



Figure 2: Tube around a two dimensional manifold. The manifold is shown in red, 
and the tube is divided into a main part, half-cylinders around the edges, and corner 
wedges. 



A version of the tube formula, without boundary corrections, was first 
derived by Wcyl (1939). Naimau (1990) provided boundary corrections. The 
result is a series with d+1 terms. The first four terms (for manifolds on the 
unit sphere) are 

Volume = ^0^"- P{B(d+iy2,(n-d-i)/2 > ™^) 

Ad+l 

+ ^JjP[Bd/2,{n^d)/2 > W ) 

K2+ll+mo An ^ 2x 

H -I -P(i!(d-l)/2,(n-d+l)/2 > W ) 

2n Ad-i 



h + mi+no An , ^ 2^ 

^ -^—P(B(d-2)/2,(n-d+2)/2 > W ) 



The constants Zo, and I2 arise from the corresponding series for the half- 
tubes around the boundaries of the manifold. is the (d — l)-dimensional 
volume of the boundaries (or the total length of the four edges in Figure 2). 
li and I2 are higher order terms representing boundary curvature. 

mo and mi arise from the 'corner wedges' where two boundary faces meet. 
In figure 2, mo is the sum of the four wedge angles at each corner of the 
manifold. 

no arises for manifolds with d > 3, from the corners where three (or more) 
boundary faces meet. 

2.1 Random Processes 

In statistical applications, the fundamental use of the tube formula is to find 
(or at least approximate) the distribution of the maximum of certain random 
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processes. As a simple example, consider the process 

Z{\) ^ {T{X),U) 

where T(A) is an /?"-valued vector function, and U is uniformly distributed 
over the unit sphere. Suppose that one is interested in finding 

P(supZ(A) > w) 

for some w. 

The inner product exceeds w if, and only if, U is sufficiently close to T{X). 
Specifically, 

ItT(A) - Uf = ||r(A)f - 2 (T(A), U) + \\Uf = 2(1 - (T(A), [/)). 

Hence, ||T(A) - U\\ < r if, and only if, (T(A), U) > w, where = 2(1 - w). 
The probability (2.1) is therefore simply 

Area of tube of radius r around {r(A)} 



Surface area of unit sphere in 7?" 

In many statistical applications, one is interested in the distribution of the 
maximum of a Gaussian process, 

where e follows the standard multivariate normal distribution. To reduce this 
to the uniform process, one needs to condition on the length of the e vector, 
and integrate over the conditional distribution; see Sun and Loader (1994). 
The final result, up to fourth order, is 

P(supZ(A)>c) ^ ^P{^f^^^^^^>c^) 

Ad+l 

h + mi + no ,2 ^ 2x 

where Xk denotes a chi-square random variable with k degrees of freedom. 



3 The libtube Library 

The main computational problem in implementing results based on the tube 
formula is evaluation of the constants kq, K2, lo e.t.c. The libtube library 
implements the tube library up to fourth order terms. To use the library, one 
must first write a 'manifold function' defining the problem, libtube takes 
the manifold function as input, and uses numerical integration methods to 
compute the constants. 

The library can be downloaded from http : //www . herine . net/stat/libtube . 
The library is written in C, and can be compiled on Linux systems using 

y, make 

y, make install 

to install the libraries in /usr/local/lib. 

y, cc -o nlreg nlreg.c -Itube -Imut -Im 

The library (and the examples given in this paper) have been written and 
tested using the Gnu C compiler available in most Linux distributions. The C 
code should be compatible with most other compilers and operating systems. 
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3.1 Manifold Functions 



Suppose the manifold is defined by a vector function T{x) mapping a d- 
dimensional domain X to the manifold A4 in n-dimensional space. The con- 
stants in the tube formula can be computed from T{x) and its derivatives, so 
in it's simplest form, the manifold function simply computes these. In sta- 
tistical applications, one usually doesn't get T{x) naturally, but rather one 
gets a vector l{x) such that T{x) = i(a::)/||Z(a;) || (see the regression examples 
in Sections 4 and 5). The manifold function can instead provide l{x) and its 
derivatives. 

In still other examples, one doesn't even obtain l{x) directly, but instead 
obtains a covariance function cr{x, x') — {l{x), l{x')) (see the mixture example. 
Section 6). Since the distance between any two points on the manifold is given 
by 

\\l{x) — l{x )\\ = a{x,x) + a{x ,x ) — 2a{x,x ), 

knowledge of the covariance function determines l{x) up to an orthogonal 
transformation. The manifold function can provide a{x, x') and its derivatives. 

The precise form of the manifold functions is illustrated by the examples. 
After writing the manifold function, the most useful functions in libtube are: 

• tube_contstants(), to numerically evaluate and the other constants 
appearing in (??). 

• tailpO and critvalO, which compute tail probabilities corresponding 
to a specified cut-off, and critical values corresponding to a specified 
significance level. 

Calling sequence for tube_constants () . 

The function to compute the constants is 

int tube_constants (f , d, n, ev, mg, fl, kap, wk, deb, uc) ; 

int (*f)(); 

int d, n, ev, mg; 

double *fl, *kap, *wk; 

int deb , uc ; 

The arguments to this function are: 

• f The manifold function to compute l(x) and its derivatives. 

• d The dimension of the manifold. 

• m The maximum length of the l{x) vectors. The argument provided is 
only used to allocate work space; the actual length of l[x) is returned by 
the manifold function. 

• ev Integration type. For rectangular domains, ISIMPSON is the most 
useful. 

• mg Integer vector, giving the number of partitions to use in each dimen- 
sion of the numerical integration rules. 

• f 1 Integration limits. A numeric vector with length 2d. The first d com- 
ponents give lower limits for each variable; the remaining d components 
give upper limits. 

• kap is the vector through which the computed constants are returned. 
It should be allocated with at least min(ci -f 1,4) terms. The values 
returned are kq, /o/2, (k2 + li + ma)/{2n) and {I2 + mi -f no)/(47r). 

• wk is a workspace vector. If wk=NULL, the required workspace will be allo- 
cated and freed within the tube_constants() function. To pre-allocate 
the space, the required length can be found by calling kO_reqd(d,m) . 
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• terms Number of terms to compute, from 1 to 4. 

• uc An indicator variable indicating whether the manifold function com- 
putes the weight vectors uc=0 or covariance derivative matrix uc=l. 

Calling sequence for tailpO and critvalO. 

Tail probabilities are computed using the function tailp: 

double tailp (c ,kO,m,d, s ,n, process) 
double c , *kO , n ; 
int m, d, s, process; 

Critical values corresponding to a specified tail probability are computed 
using the critval function: 

double cr itval (alpha, kO ,m,d, s ,n, process) 
double *kO, al, n; 
int m, d, it, s; 

These arguments represent: 

• c Cut-off value for tailp (). 

• alpha tail probability for critvalO. 

• kO is the vector of constants computed by the tube_constaiits() func- 
tion. 

• m is the number of terms in kO. This is the returned value of tube_constants () , 
and is equal to min(ci -I- 1,4). 

• d is the dimension of the manifold. 

• c critical value (tailp () only). 

• s Either ONE_SIDED or TWO_SIDED. 

• n For the t-process, the residual degrees of freedom used to estimate 
a. For the uniform process, the dimension n. Ignored for the Gaussian 
process. Beware that n must have type double. 

• process Either GAUSS (when e is multivariate Gaussian); TPROC (Gaus- 
sian process with estimated variance); or UNIF (when e is uniform on the 
unit sphere). 

3.2 Writing a Manifold Function with vectors 

The manifold function computes the vector l{x) and its derivatives. The basic 
form of the function is 

int mymf (x, 1 ,reqd) 
double *x, *1; 
int reqd; 

-[ /* function body goes here */ 
> 

The x argument is a point in the input space; 1 is a vector to be filled in by 
the manifold function. The final argument, reqd, is an integer indicating what 
the library requires from the manifold function. If reqd=0, only the vector 
l{x) is required. If reqd=l, then both l{x) and l'{x) (or all the first-order 
partial derivative vectors of l{x)) are required. If reqd=2, then additionally 
the second-order partial derivative vectors are required. 

In statistical applications, the manifold function will generally require a 
data vector, sample size n, dimension d and variables other than x in order to 
perform its calculations. These variables should be assigned to global variables 
so that they are accessible in the manifold function. 
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The results of the computations are returned through the 1 vector. The 
vector l{x) is placed in the first n elements. The first-order derivatives are 
placed in the next n x d elements. The second-order derivatives are placed in 
the next n x d x d elements. 

The function should return n, the length of the vector l{x). Generally, 
this should be equal to the n value provided in the tube_constants() call; 
it should never be larger. It can be less. An example where it may be less 
is for a kernel regression with compactly supported kernel; only the non-zero 
elements of l{x) need be retained. 



3.3 Writing a Manifold Function with a covariance 
function. 

The structure of a manifold function based on the covariance is identical to 
the vector case; it differs in what is computed. 

Given a covariance function a{x,x'), the manifold function needs to com- 
pute (in the one-dimensional case), 

/ /\ d(T(x,x') d'^a(x,x') 

g{x,x) — aL/^ ' 

d(T{x ,x' ) (t{xjx') a{x ,x' ) 

dx dxdx' dxdx'"^ 

d'^ ^{x ,x' ) ^{xjx') d'^<j{x,x' ) 

dx'-^ dx-^dx' dx'-^dx''^ . 

evaluated at = x. Again, the matrix is stored in the vector 1, with the 
columns stacked atop each other. 

In higher dimensions, the required matrix is most easily written in terms 
of differential operators. The required {1 -\- d -\- d^) x + + matrix is 

/ ^ \ 

Dx, 



Dx 
Dx, 



a{x,x') (I D^ 



^ x'^ ^ x'^,x'^ 



D^ 



\ Dx^,x^ / 

where D represents the partial derivative operator with respect to the sub- 
scripted variables. 

Another view is as follows. If L is the matrix computed by a manifold func- 
tion with vectors, then L"^L is the matrix computed by a manifold function 
with a covariance function. 



4 Example: Testing in Nonlinear Regres- 
sion 

This was the motivating example for Hotclhng (1939), and was developed in 
much more detail by Knowles and Sicgmund (1989). Suppose one has data 
{xi,Yi),i = 1, . . . , n, and a nonlinear regression model, such as 

Y,^ae^^'+e,. (1) 

The important feature of this model is that the parameter a enters the model 
linearly, while 7 enters nonlinearly. Assume that the errors are independent 
N{0,a^). 
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Consider the problem of testing Ho : a = vs Hi : a ^ 0. It can be shown 
that the log-likelihood ratio test statistic is equivalent to 



_ min^,^||y-a/(7)|j 

where = (e^"i , . . . , e^"" ). 

In classical statistical theory, log-likeliood ratio statistics often have asymp- 
totic distributions. However, this is not the case for the statistic (2). One 
way to see this is to recall that proofs of the results are based on a quadratic 
expansion of the statistic under the null parameters. For the present problem 
this would require an expansion around (0, 70) where 70 is 'the' null value of 
7. Unfortunately this is undefined: when q = 0, the parameter 7 does not 
appear in (1); it is not identifiable! 

For fixed 7, minimizing over a is a linear least-sqaures problem. It follows 
that 

■ lil) Y 



L — 1 — sup 



ii'(7)ii' ri 



The null hypothesis Ho is rejected if L < 1—w for some w > 0, or equivalently, 
if 

K7) Y 



sup 



> w. 



\m\\'\\Y\\ 

The constant w must be chosen to obtain a specified significance level. That 
is, we need to be able to evaluate probabilities of the form 

P(sup|(r(7),f/)|>u,) (3) 

7 

where T{'y) = ^(7)/||K7)ll defines a curve on the unit sphere, and U = y/||y|| 
is (under Ho : a = 0) uniformly distributed on the surface of the sphere. 



4.1 Non-linear Regression: Implementation 

Code implementing the tube formla for the non-linear regression problem is 
shown below. The program consists of the manifold function regmf () , and the 
main routine mainO that reads in the data and computes the tube constants. 
Note that the data vectors and sample size are stored as global variables, so 
that they can be accessed within the manifold function. 

The manifold function computes the components of ^(7); li{'y) = e^^', and 
of ^'(7), I'ii'y) = "fe'^^' . These vectors are stored end-to-end in the 1 argument. 

#include <stdio.h> 
#include <math . h> 
#include <tube . h> 
#define MAXN 1000 



double x[MAXN], y[MAXN] ; 
int n; 

int regmf (gam, 1 ,reqd) 
double *gain, *1; 
int reqd; 
■[ int i ; 

double *11 ; 

11 = &l[n] ; 

for (i=0; i<n; i++) 

{ l[i] = exp (gam [0] *x [i] ) ; 
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11 [i] = x[i]*exp(gain[0]*x[i]) ; 

y 

return (n) ; 

} 

int mainO 

{ FILE *infile; 

char filename [100] ; 

int i , mg ; 

double gamlimits [2] , kappa [4] ; 

printfC'Data filename ? "); scanf (""/s" , filename) ; 
printfC'n = ? ") ; scanf (""/.d" ,&n) ; 
infile = f open(f ilename, "r") ; 

for (i=0; i<n; i++) fscanf (infile , ""/.If '/.If " ,&x [i] ,&y [i] ) ; 
gamlimits [0] = -2.0; 
gamlimits [1] = 2.0; 
mg = 100; 

tube_constants (regmf , 1 , n , ISIMPSON , &mg , gamlimits , kappa , NULL ,0,0); 
printf ("y.8.5f 7.8. 5f \n" , kappa [0] ,kappa[l]); 

} 



5 Example: Simultaneous Confidence Bands 

Application of the tube formula to find simultaneous confidence bands for 
regression models has been studied in Naiman (1987), Sun and Loader (1994) 
among others. Consider again regression data, but now suppose that the 
model is 

Yi = ao + aiXi + a2xl + e, = /i(a;i) + 

(although we formulate the problem for quadratic regression, extension to 
other linear models is straightforward). The goal is to find confidence bands 

jl{x) ± c\/ var(/i(x)) 
with simultaneous coverage over some nice domain X: 

P{\fi{x) - ^t(x)| < ca\\l{x)\\ for all x G A:") = 1 - a. (4) 
The least-squares estimates of the parameters are 



= (X^X)"^X'^y 



where X is the design matrix. For fixed x, ij.{x) is estimated by 

fl{x) = ao + aix + a2x'^ ^ {1 x x'^ ) {X.'^Xy^X'^Y ^ {l(x),Y) , 

and the variance of the estimate is var(/i(a;)) — a'^\\l{x)\\. 

Now, fi{x) — nix) = {l{x),e}, and the probability (4) is equivalent to 

This problem can be solved using the Gaussian process variant of the tube 
problem. 
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Suppose X = QR is the Q^-decomposition of the design matrix. Then 
l{x) lies in the column space of Q for all x, and so 

so it suffices to work with the vector I'ix) — Cl^l{x) = (R-^)^^/(a::) where 
f{x) is a vector of the polynomial basis functions. The derivatives are easily 
found; 



Ar(.).(R^)-A/(.) 



and so on. 



5.1 Simultaneous Confidence Bands: Implementa- 
tion 

Code for the quadratic regression computations, in an arbitrary number of 
dimensions, is shown below. The functions quadO, quadiO and quadij () 
compute the quadratic basis functions f{x); first-order partial derivatives 
and second-order partial derivatives respectively. The manifold function is 
quadmf (). The mainO function reads the data from a file; computes the de- 
sign matrix and its QR-decomposition; and then calls the tube_constants() 
function (The QR functions, qr() and qrtinvxO, as well as transpose (), 
are part of the mut library). 

When the program is run, the user is prompted for a data file (containing 
a matrix of the predictor variables); data dimension (n and d), and limits for 
the confidence band computation. 

The tube constants are computed, then the critical value c for 95% confi- 
dence bands. Note that the final argument to critval is the residual degrees 
of freedom used to estimate a; Sun and Loader (1994) give the modification 
of (??) for this case. 

#include <stdio.h> 
#include <tube . h> 
#include <mutil.h> 



int dim, n, p; 
double *X ; 



void quad(x,f) 
double *x, *f ; 
{ int i , j , k ; 
k = 0; 

f[k++] = 1.0; 

for (i=0; Kdim; i++) f [k++] = x[i]; 
for (i=0; Kdim; i++) 
for (j=i; j<dim; j++) 
f [k++] = x[i]*x[j] ; 

} 



void quadi (x,f , iO) 
double *x, *f ; 
int iO; 

-[ int i , j , k ; 
k = 0; 

f[k++] = 0.0; 
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for (i=0; i<dim; i++) f [k++] = (i==iO) ; 
for (i=0; i<dim; i++) 
for (j=i; j<dim; j++) 

f[k++] = (i==iO)*x[j] + (j==iO)*x[i] ; 

} 

void quadij (x,f ,iO, jO) 
double *x, *f ; 
int iO, jO; 
{ int i , j , k ; 
k = 0; 

f[k++] = 0.0; 

for (i=0; i<dim; i++) f [k++] =0.0; 
for (i=0; i<dim; i++) 
for (j=i; j<dim; j++) 

f[k++] = ((i==iO) & (j==jO)) + ((i==jO) & (j==iO)); 

> 

int quadmf (x , 1 , reqd) 

double *x, *1; 

int reqd; 

{ int i , j , k ; 

k = 0; 

quad(x,l) ; 

qrtinvx(X,l,n,p) ; 

k++; 

for (i=0; i<dim; i++) 
{ quadi(x,&l [k*p] ,i) ; 

qrtinvx(X,&l [k*p] ,n,p) ; 

k++; 

} 

for (i=0; i<dim; i++) 
for (j=0; j<dim; j++) 
{ quadij (x,&l [k*p] ,i,j) ; 
qrtinvx(X,&l [k*p] ,n,p) ; 
k++; 

y 

return (p) ; 

> 

int mainO 

{ FILE *infile; 

char filename [100] ; 

int i, j , mg[100] ; 

double xlim[100], kappa [4] , datarow [100] ; 
printfC'Data filename ? "); scanf ("7,3" , filename) ; 
printf("n = ? "); scanf ("7.d" ,&n) ; 
printf("dim = ? "); scanf ("'/.d" ,Mim) ; 
infile = fopen(f ilename, "r") ; 
if (inf ile==NULL) 

{ printf ("Error : can't read input file\n"); 
return (0) ; 

} 

p = 1 + dim + dim*(dim+l)/2; 

X = (double *) calloc (n*p, sizeof (double) ) ; 
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for (i=0; i<n; i++) 

{ for (j=0; j<dim; j++) 

f scanf ( inf ile , " "/.If " , &dat arow [ j ] ) ; 
quad(datarow,&X[i*p] ) ; 

y 

transpose (X,n,p) ; 
qr(X,n,p,NULL) ; 

for (i=0; i<dim; i++) mg[i] = 20; 
xlim[0] = -2; xlim[l] = -2; 
xlim[2] = 2; xlim[3] = 2; 

tube_constants (quadmf , dim, p, IS IMPSON,mg,xlim, kappa, NULL, 3,0) ; 

printf ("kappa: 7.8. 5f 7.8. 5f 7.8 . 5f 7.8. 5f \n" , kappa [0] , kappa [1] , kappa [2] , kappa [3] ) ; 



6 Example: Mixture Models 

Suppose Xi , ■ ■ ■ , Xn are an i.i.d. sample from a density 

fa,\{x) = (1 - a)fo{x) + a(l>{x,X) 

where a and A are unknown parameters, with < a < 1. The object is to 
test Ho : a = vs Hi : a > 0. This is a simple example of mixture testing: 
under Ho, the single component /o(a:) describes the data, while under Hi, the 
two components are required. 

Consider the normalized score process proposed by Pilla and Loader (2003). 
The score process is 



5(A) = ■ 



z — 1 

Under the null hypothesis, this has mean and covariance function na{\, A^), 
where 

,(A,At)= / 



The normalized score process is S*(A) = S(A) / na{X, A). This asymptot- 
ically behaves like a Gaussian process Z{X), with mean under Ho, and 
a nonzero mean under Hi. The maximum of the normalized score process 
serves as the test statistic. 

Since an explicit vector representation of Z{\) is not readily available, 
the manifold function (mixmf in the code below) must be written using the 
covariance function and its partial derivatives. 

There is one additional difficulty. The normalized score process has a 
singularity at = 0. For this reason, the manifold function works with Taylor 
series expansions of the covariance in this reason. Also, the singularity results 
in a discontinuity in S*{X), and lo — 4. 

The main routine in the program below sets limits for /i, calls the tube_constants () 
function, and computes the 5% critical value. 

kappaO = 5.27449 
10/2 = 2.00000 
Level 0.05 critical value = 2.49455 

#include <stdio.h> 
#include <math . h> 
#include <tube . h> 
#define MAXN 1000 
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double x[MAXN], y[MAXN] ; 
int n; 

int mixmf (mu, 1 ,reqd) 
double *mu, *1; 
int reqd; 
-[ double emm, mm; 

if (fabs(mu[0]) < 0.01) 
{ mm = mu [0] *mu [0] ; 

1[0] = 1 + mm/2*(l + mm/3*(l + mm/4*(l + mm/5))); 

1[1] = 0.5*(1 + mm/3* (2 + mm/4* (3 + mm/5* (4 + 5/6*mm)))); 

1[1] = 1[2] = mu[0]*l[l] ; 

1[3] = 0.5*(1 + mm/3* (4 + mm/4* (9 + mm/5* (16+25/6*mm) ) ) ) ; 
} else 

{ emm = exp (mu [0] *mu [0] ) ; 
1[0] = emm-1; 
1[1] = 1[2] = mu[0]*emm; 
1[3] = emm*(l+mu[0] *mu[0] ) ; 

} 

return (2) ; 

> 

int mainO 

■[ int i , mg , t ; 

double mulimits[2], kappa [4] ; 

mulimits[0] = -3.0; 

mulimits[l] = 3.0; 

mg = 200; 

t = tube_const ants (mixmf , 1 ,n, ISIMPSON.&mg.mulimits , kappa, NULL ,2,1); 
/* modify kappa [1] = 10/2 for the singularity */ 
kappa [1] += 1.0; 

printf ("kappaO = 7.8 . 5f \n" ,kappa[0] ) ; 
printfC 10/2 = •/.8.5f \n" ,kappa[l] ) ; 

printf ("Level 0.05 critical value = 7.8 . 5f \n" , critvaKkappa, t , 1 , . 05, 10 , 1 , . 0) ) ; 

> 

References 

Hotelling, H. (1939). Tubes and spheres in ri-spaces, and a class of statistical 
problems. American Journal of Mathematics, 61:440-460. 

Knowles, M. and Siegmund, D. (1989). On Hotelling's geometric approach 
to testing for a nonlinear parameter in regression. International Statistical 
Review, 57:205-220. 

Naiman, D. Q. (1987). Simultaneous confidence bounds in multiple regression 
using predictor variable constraints. Journal of the American Statistical 
Association, 82:214-219. 

Naiman, D. Q. (1990). On volumes of tubular neighborhoods of spherical 
polyhedra and statistical inference. The Annals of Statistics, 18:685-716. 

Pilla, R. S. and Loader, C. (2003). The volume-of-tube formula: Perturbation 
tests, mixture models and scan statistics. Unpublished Manuscript. 



13 



Sun, J. and Loader, C. (1994). Simultaneous confidence bands for linear 
regression and smoothing. The Annals of Statistics, 22:1328-1345. 

Weyl, H. (1939). On the volume of tubes. American Journal of Mathematics, 
61:461-472. 



14 



